Blame_analysis
Markus Lundsfryd Jensen and Rune Trust
2025-11-04
#load packages
#assumption check function
overdisp_fun <- function(model) {
rdf <- df.residual(model)
rp <- residuals(model,type="pearson")
Pearson.chisq <- sum(rp^2)
prat <- Pearson.chisq/rdf
pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}
model_check <- function(model_name){
hist(resid(model_name))
overdisp_fun(model_name)
plot(fitted(model_name), resid(model_name)) #residuals vs fitted
abline(h=0)
sim_res <- simulateResiduals(model_name, n = 1000)
# QQ plot
plotQQunif(sim_res) # QQ plot for simulated residuals
plot(sim_res)
testZeroInflation(sim_res)
print(overdisp_fun(model_name))
}#load data
## Rows: 2191051 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): text, speaker, party, parties_in_government
## dbl (4): paragraph, sentence_nr, prediction, confidence
## lgl (1): current_speaker_in_government
## date (1): date
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.#show daily/montly counts of blame pr party and in total #Show daily/monthly relative amounts of blame pr party #Both of the above show also in regards to the amount of sentences by the party
#Put in lines of personal political scandales and global times of crisis
#Make model showing differences in amount of blame for center parties vs right/left wing #The same for opposition do the opposition always shwo more blame (here also interaction between right/left wing) (glmer with random intercept?)
## `summarise()` has grouped output by 'month'. You can override using the
## `.groups` argument.#Monthly
ggplot(monthly_blame, aes(x = month, y = total_blame, color = party)) +
geom_smooth(span = 0.5, se = TRUE) + # per-party trends
geom_smooth(aes(x = month, y = total_blame), # global trend
color = "black", linetype = "dashed",
se = FALSE, span = 0.5) +
labs(
title = "Monthly blame Over Time by Party with Global Trend",
x = "Month",
y = "Total Blame"
) +
theme_minimal()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : span too small. fewer data values than degrees of freedom.## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : at 2022.4## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : radius 6.9444e-07## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : all data on boundary of neighborhood. make span bigger## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 2022.4## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 0.00083333## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 1## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : at 2022.6## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : radius 6.9444e-07## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : all data on boundary of neighborhood. make span bigger## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 6.9444e-07## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : zero-width neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : zero-width neighborhood. make span bigger## Warning: Failed to fit group 2.
## Caused by error in `predLoess()`:
## ! NA/NaN/Inf in foreign function call (arg 5)## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
#DO the above as a proportion of sentences uttered
#Lidstone smoothing, additive smoothing
monthly_blame <- monthly_blame %>%
mutate(offset_ratio = (total_blame+0.001)/(total_nr_sentences+0.001))ggplot(monthly_blame, aes(x = month, y = offset_ratio, color = party)) +
geom_smooth(span = 0.5, se = FALSE) + # per-party trends
geom_smooth(aes(x = month, y = offset_ratio), # global trend
color = "black", linetype = "dashed",
se = FALSE, span = 0.5) +
labs(
title = "Monthly blame Over Time by Party with Global Trend",
x = "Month",
y = "Blame pr sentence"
) +
theme_minimal()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : span too small. fewer data values than degrees of freedom.## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : at 2022.4## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : radius 6.9444e-07## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : all data on boundary of neighborhood. make span bigger## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 2022.4## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 0.00083333## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 1## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : at 2022.6## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : radius 6.9444e-07## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : all data on boundary of neighborhood. make span bigger## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 6.9444e-07## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : zero-width neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : zero-width neighborhood. make span bigger## Warning: Failed to fit group 2.
## Caused by error in `predLoess()`:
## ! NA/NaN/Inf in foreign function call (arg 5)## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
#Show number of sentences pr year
ggplot(monthly_blame, aes(x = month, y = total_nr_sentences))+
geom_smooth(aes(x = month, y = total_nr_sentences), span = 0.5, se = FALSE)+
labs(
title = "Number of uttered sentences monthly",
x = "Month",
y = "Total sentences"
) +
theme_minimal()## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
In order to make it easier to do statistics and interpretation.Focus will not be on specific parties, but on left/right and center/wing
right_wing_parties <- c("V", "KF", "LA", "KD", "DF", "NB", "DD")
center_parties <- c("S", "V", "KF")
monthly_blame$block <- ifelse(monthly_blame$party %in% right_wing_parties, "Right", "Left")
monthly_blame$wing <- ifelse(monthly_blame$party %in% center_parties, "Center", "Wing")
monthly_blame$wing <- as.factor(monthly_blame$wing)
monthly_blame$block <- as.factor(monthly_blame$block)
monthly_blame$in_gov <- as.factor(monthly_blame$in_gov)
monthly_blame$month_index <- 12 * (as.numeric(format(monthly_blame$month, "%Y")) - 2000) + (as.numeric(format(monthly_blame$month, "%m")) - 1)str(monthly_blame) #all good## tibble [2,330 × 10] (S3: tbl_df/tbl/data.frame)
## $ month : 'yearmon' num [1:2330] Jan 2000 Jan 2000 Jan 2000 Jan 2000 ...
## $ party : Factor w/ 30 levels "-","ALT","CD",..: 5 6 13 14 22 23 24 30 5 6 ...
## $ in_gov : Factor w/ 2 levels "FALSE","TRUE": 1 1 1 1 2 2 1 1 1 1 ...
## $ total_nr_sentences: int [1:2330] 817 585 234 259 815 993 469 303 715 783 ...
## $ total_blame : num [1:2330] 67 26 16 20 40 41 34 15 18 57 ...
## $ year : num [1:2330] 0 0 0 0 0 0 0 0 0 0 ...
## $ offset_ratio : num [1:2330] 0.082 0.0444 0.0684 0.0772 0.0491 ...
## $ block : Factor w/ 2 levels "Left","Right": 2 1 2 2 1 1 1 2 2 1 ...
## $ wing : Factor w/ 2 levels "Center","Wing": 2 2 2 1 2 1 2 1 2 2 ...
## $ month_index : num [1:2330] 0 0 0 0 0 0 0 0 1 1 ...#Blame predicted by time
model_blame_time_0 <- glmmTMB(total_blame ~ month_index + (1|in_gov) + offset(log(total_nr_sentences)),
ziformula = ~1,
family = nbinom2,
data = monthly_blame)model_blame_time_1 <- glmmTMB(total_blame ~ month_index + (1|party) + offset(log(total_nr_sentences)),
ziformula = ~1,
family = nbinom2,
data = monthly_blame)model_blame_time_2 <- glmmTMB(total_blame ~ month_index + (1|party) + (1|in_gov) + offset(log(total_nr_sentences)),
ziformula = ~1,
family = nbinom1,
data = monthly_blame)
summary(model_blame_time_2)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ month_index + (1 | party) + (1 | in_gov) + offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 18945.4 18979.9 -9466.7 18933.4 2324
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.07621 0.2761
## in_gov (Intercept) 0.04032 0.2008
## Number of obs: 2330, groups: party, 13; in_gov, 2
##
## Dispersion parameter for nbinom1 family (): 5.54
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.7610684 0.1646103 -16.773 < 2e-16 ***
## month_index -0.0008101 0.0000991 -8.175 2.97e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.67 1151.57 -0.019 0.985anova(model_blame_time_0, model_blame_time_1, model_blame_time_2)## Data: monthly_blame
## Models:
## model_blame_time_0: total_blame ~ month_index + (1 | in_gov) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_blame_time_1: total_blame ~ month_index + (1 | party) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_blame_time_2: total_blame ~ month_index + (1 | party) + (1 | in_gov) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## model_blame_time_0 5 19254 19283 -9622.2 19244
## model_blame_time_1 5 19133 19162 -9561.4 19123 121.49 0 < 2.2e-16 ***
## model_blame_time_2 6 18945 18980 -9466.7 18933 189.54 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Plot best model
#install.packages("sjPlot")
library(sjPlot)
sjPlot::plot_model(model_blame_time_2, type = "re")## [[1]]
##
## [[2]]
#Due to the offset, in order to do meaningfull interpretation, the total number of sentences must be kept constant. Thus the following plot shows the change in predicted amount of blame pr 1000 sentences uttered as time goes on
sjPlot::plot_model(model_blame_time_2,
type = "pred",
terms = "month_index",
condition = c(total_nr_sentences = 100),
axis.title = c("Months after Jan 2000", "Predicted Blame Counts pr 100 Uttered Sentences"),
title = "General Evolution of Blame in the Danish Parliament over time",
ci.lvl = NA,
show.legend = TRUE
) +
theme_minimal()## You are calculating adjusted predictions on the population-level (i.e.
## `type = "fixed"`) for a *generalized* linear mixed model.
## This may produce biased estimates due to Jensen's inequality. Consider
## setting `bias_correction = TRUE` to correct for this bias.
## See also the documentation of the `bias_correction` argument.
#shows downtrend over time
#Create models of overall tendentcy since 2000 of increasing complexity and anova them
model_intercept <- glmmTMB(total_blame ~ 1 + (1|party) + (1|month_index)+ offset(log(total_nr_sentences)),
ziformula = ~1,
family = nbinom1,
data = monthly_blame)
summary(model_intercept)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ 1 + (1 | party) + (1 | month_index) + offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 18734.1 18762.8 -9362.0 18724.1 2325
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.07830 0.2798
## month_index (Intercept) 0.05223 0.2285
## Number of obs: 2330, groups: party, 13; month_index, 274
##
## Dispersion parameter for nbinom1 family (): 3.8
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.77855 0.08258 -33.65 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.08 955.59 -0.022 0.982model_in_gov <- glmmTMB(total_blame ~ in_gov + (1|party) + (1|month_index)+ offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = monthly_blame)
summary(model_in_gov)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + (1 | party) + (1 | month_index) + offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 18410.4 18444.9 -9199.2 18398.4 2324
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.06295 0.2509
## month_index (Intercept) 0.05298 0.2302
## Number of obs: 2330, groups: party, 13; month_index, 274
##
## Dispersion parameter for nbinom1 family (): 3.14
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.74565 0.07445 -36.88 <2e-16 ***
## in_govTRUE -0.37098 0.01985 -18.69 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.51 1038.23 -0.021 0.983model_block <- glmmTMB(total_blame ~ in_gov + block + (1|party) + (1|month_index) + offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = monthly_blame)
summary(model_block)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + block + (1 | party) + (1 | month_index) +
## offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 18412.4 18452.6 -9199.2 18398.4 2323
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.06295 0.2509
## month_index (Intercept) 0.05298 0.2302
## Number of obs: 2330, groups: party, 13; month_index, 274
##
## Dispersion parameter for nbinom1 family (): 3.14
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.75539 0.10680 -25.800 <2e-16 ***
## in_govTRUE -0.37098 0.01985 -18.688 <2e-16 ***
## blockRight 0.01850 0.14572 0.127 0.899
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.66 1117.05 -0.019 0.985model_wing <- glmmTMB(total_blame ~ in_gov + wing + (1|party) + (1|month_index)+ offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = monthly_blame)
summary(model_wing)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + wing + (1 | party) + (1 | month_index) +
## offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 18410.7 18451.0 -9198.3 18396.7 2323
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.05522 0.2350
## month_index (Intercept) 0.05300 0.2302
## Number of obs: 2330, groups: party, 13; month_index, 274
##
## Dispersion parameter for nbinom1 family (): 3.14
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.90427 0.13696 -21.20 <2e-16 ***
## in_govTRUE -0.37002 0.01987 -18.62 <2e-16 ***
## wingWing 0.20962 0.15764 1.33 0.184
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.28 922.40 -0.023 0.982model_wing_block <- glmmTMB(total_blame ~ in_gov + wing + block + (1|party) + (1|month_index) + offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = monthly_blame)
summary(model_wing_block)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + wing + block + (1 | party) + (1 | month_index) +
## offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 18412.6 18458.6 -9198.3 18396.6 2322
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.05484 0.2342
## month_index (Intercept) 0.05300 0.2302
## Number of obs: 2330, groups: party, 13; month_index, 274
##
## Dispersion parameter for nbinom1 family (): 3.14
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.93644 0.16468 -17.831 <2e-16 ***
## in_govTRUE -0.36996 0.01987 -18.619 <2e-16 ***
## wingWing 0.21853 0.15923 1.372 0.170
## blockRight 0.04827 0.13826 0.349 0.727
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.38 974.08 -0.022 0.982model_wing_block_int <- glmmTMB(total_blame ~ in_gov + wing + block + wing*block + (1|party) + (1|month_index) + offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = monthly_blame)
summary(model_wing_block_int)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + wing + block + wing * block + (1 | party) +
## (1 | month_index) + offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 18414.0 18465.8 -9198.0 18396.0 2321
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.05221 0.2285
## month_index (Intercept) 0.05299 0.2302
## Number of obs: 2330, groups: party, 13; month_index, 274
##
## Dispersion parameter for nbinom1 family (): 3.14
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.81219 0.22947 -12.255 <2e-16 ***
## in_govTRUE -0.36992 0.01987 -18.618 <2e-16 ***
## wingWing 0.06676 0.25273 0.264 0.792
## blockRight -0.13827 0.28049 -0.493 0.622
## wingWing:blockRight 0.24231 0.32004 0.757 0.449
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.64 1105.78 -0.02 0.984model_wing_block_int_2 <- glmmTMB(total_blame ~ in_gov + block + in_gov*block + (1|party) + (1|month_index) + offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = monthly_blame)
summary(model_wing_block_int_2)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + block + in_gov * block + (1 | party) +
## (1 | month_index) + offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 18401.7 18447.7 -9192.8 18385.7 2322
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.06153 0.2481
## month_index (Intercept) 0.05271 0.2296
## Number of obs: 2330, groups: party, 13; month_index, 274
##
## Dispersion parameter for nbinom1 family (): 3.12
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.748347 0.105649 -26.014 < 2e-16 ***
## in_govTRUE -0.444119 0.028757 -15.444 < 2e-16 ***
## blockRight -0.001707 0.144255 -0.012 0.990559
## in_govTRUE:blockRight 0.177978 0.050112 3.552 0.000383 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.37 956.14 -0.022 0.982model_wing_block_int_3 <- glmmTMB(total_blame ~ in_gov + wing + block + in_gov*block + (1|party) + (1|month_index) + offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = monthly_blame)
summary(model_wing_block_int_3)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + wing + block + in_gov * block + (1 | party) +
## (1 | month_index) + offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 18401.7 18453.5 -9191.9 18383.7 2321
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.05280 0.2298
## month_index (Intercept) 0.05272 0.2296
## Number of obs: 2330, groups: party, 13; month_index, 274
##
## Dispersion parameter for nbinom1 family (): 3.12
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.93526 0.16169 -18.154 < 2e-16 ***
## in_govTRUE -0.44354 0.02876 -15.423 < 2e-16 ***
## wingWing 0.22576 0.15634 1.444 0.148715
## blockRight 0.02874 0.13591 0.211 0.832531
## in_govTRUE:blockRight 0.17926 0.05013 3.576 0.000349 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.54 1041.90 -0.021 0.984model_wing_block_int_4 <- glmmTMB(total_blame ~ in_gov + wing + block + wing*block + in_gov*wing + in_gov*block + (1|party) + (1|month_index) + offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = monthly_blame)
summary(model_wing_block_int_4)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + wing + block + wing * block + in_gov *
## wing + in_gov * block + (1 | party) + (1 | month_index) +
## offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 18404.6 18467.9 -9191.3 18382.6 2319
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.04898 0.2213
## month_index (Intercept) 0.05271 0.2296
## Number of obs: 2330, groups: party, 13; month_index, 274
##
## Dispersion parameter for nbinom1 family (): 3.12
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.78716 0.22250 -12.526 < 2e-16 ***
## in_govTRUE -0.43716 0.03276 -13.344 < 2e-16 ***
## wingWing 0.04535 0.24511 0.185 0.853202
## blockRight -0.19506 0.27226 -0.716 0.473716
## wingWing:blockRight 0.29039 0.31057 0.935 0.349779
## in_govTRUE:wingWing -0.02709 0.05868 -0.462 0.644287
## in_govTRUE:blockRight 0.17416 0.05273 3.303 0.000956 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.55 1047.29 -0.021 0.984#Check assumptions for all models
model_check(model_block)
## chisq ratio rdf p
## 2083.921503 0.897082 2323.000000 0.999856model_check(model_wing)
## chisq ratio rdf p
## 2082.8294582 0.8966119 2323.0000000 0.9998653model_check(model_wing_block)
## chisq ratio rdf p
## 2082.8438645 0.8970042 2322.0000000 0.9998571model_check(model_wing_block_int)
## chisq ratio rdf p
## 2083.0641701 0.8974856 2321.0000000 0.9998466anova(model_intercept, model_in_gov, model_wing, model_wing_block, model_wing_block_int,model_wing_block_int_2,model_wing_block_int_3,model_wing_block_int_4)## Data: monthly_blame
## Models:
## model_intercept: total_blame ~ 1 + (1 | party) + (1 | month_index) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_in_gov: total_blame ~ in_gov + (1 | party) + (1 | month_index) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_wing: total_blame ~ in_gov + wing + (1 | party) + (1 | month_index) + , zi=~1, disp=~1
## model_wing: offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_wing_block: total_blame ~ in_gov + wing + block + (1 | party) + (1 | month_index) + , zi=~1, disp=~1
## model_wing_block: offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_wing_block_int_2: total_blame ~ in_gov + block + in_gov * block + (1 | party) + , zi=~1, disp=~1
## model_wing_block_int_2: (1 | month_index) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_wing_block_int: total_blame ~ in_gov + wing + block + wing * block + (1 | party) + , zi=~1, disp=~1
## model_wing_block_int: (1 | month_index) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_wing_block_int_3: total_blame ~ in_gov + wing + block + in_gov * block + (1 | party) + , zi=~1, disp=~1
## model_wing_block_int_3: (1 | month_index) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_wing_block_int_4: total_blame ~ in_gov + wing + block + wing * block + in_gov * , zi=~1, disp=~1
## model_wing_block_int_4: wing + in_gov * block + (1 | party) + (1 | month_index) + , zi=~1, disp=~1
## model_wing_block_int_4: offset(log(total_nr_sentences)), zi=~1, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df
## model_intercept 5 18734 18763 -9362.0 18724
## model_in_gov 6 18410 18445 -9199.2 18398 325.6978 1
## model_wing 7 18411 18451 -9198.3 18397 1.6774 1
## model_wing_block 8 18413 18459 -9198.3 18397 0.1217 1
## model_wing_block_int_2 8 18402 18448 -9192.8 18386 10.8854 0
## model_wing_block_int 9 18414 18466 -9198.0 18396 0.0000 1
## model_wing_block_int_3 9 18402 18454 -9191.9 18384 12.2852 0
## model_wing_block_int_4 11 18405 18468 -9191.3 18383 1.0781 2
## Pr(>Chisq)
## model_intercept
## model_in_gov <2e-16 ***
## model_wing 0.1953
## model_wing_block 0.7272
## model_wing_block_int_2 <2e-16 ***
## model_wing_block_int 1.0000
## model_wing_block_int_3 <2e-16 ***
## model_wing_block_int_4 0.5833
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(model_intercept, model_in_gov,model_wing_block_int_2,model_wing_block_int_3)## Data: monthly_blame
## Models:
## model_intercept: total_blame ~ 1 + (1 | party) + (1 | month_index) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_in_gov: total_blame ~ in_gov + (1 | party) + (1 | month_index) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_wing_block_int_2: total_blame ~ in_gov + block + in_gov * block + (1 | party) + , zi=~1, disp=~1
## model_wing_block_int_2: (1 | month_index) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_wing_block_int_3: total_blame ~ in_gov + wing + block + in_gov * block + (1 | party) + , zi=~1, disp=~1
## model_wing_block_int_3: (1 | month_index) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df
## model_intercept 5 18734 18763 -9362.0 18724
## model_in_gov 6 18410 18445 -9199.2 18398 325.6978 1
## model_wing_block_int_2 8 18402 18448 -9192.8 18386 12.6845 2
## model_wing_block_int_3 9 18402 18454 -9191.9 18384 1.9602 1
## Pr(>Chisq)
## model_intercept
## model_in_gov < 2e-16 ***
## model_wing_block_int_2 0.00176 **
## model_wing_block_int_3 0.16150
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#check assumptions of best model
model_check(model_wing_block_int_2)
## chisq ratio rdf p
## 2081.4510575 0.8964044 2322.0000000 0.9998688#SOme deviation, but cannot be made better
#plot best model (Intercept) -2.748347 0.105649 -26.014 < 2e-16
in_govTRUE -0.444119 0.028757 -15.444 < 2e-16
blockRight -0.001707 0.144255 -0.012 0.990559
in_govTRUE:blockRight 0.177978 0.050112 3.552 0.000383 ***
sjPlot::plot_model(model_wing_block_int_2,
type = "int",
#terms = c("in_gov", "block"),
condition = c(total_nr_sentences = 100),
axis.title = c("In government", "Predicted Blame Counts pr 100 Uttered Sentences"),
title = "Blame Predicted by Government and Block",
ci.lvl = NA,
show.legend = TRUE,
axis.lim = c(3.5,7)
) +
theme_minimal()## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
sjPlot::plot_model(model_wing_block_int_2,
type = "pred",
terms = c("block"),
condition = c(total_nr_sentences = 100),
axis.title = c("In government", "Predicted Blame Counts pr 100 Uttered Sentences"),
title = "Blame Predicted by Block",
ci.lvl = NA,
axis.lim = c(3.5,7)
) +
theme_minimal()## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
#The first plot shows clear differences in government along with the positive interaction term as illustrated by the different nature of black by the levels of government
#The second plot shows the lack of influence by block alone as found by the model
#Something seems to happen (see the time_series plot) in recent years, thus data since 2019 will also be evaluated (nye Borgerlige enters the room)
#Make data
monthly_blame_18 <- monthly_blame %>%
filter(year >=19)#plot data
ggplot(monthly_blame_18, aes(x = month, y = offset_ratio, color = party)) +
geom_smooth(span = 0.5, se = FALSE) + # per-party trends
geom_smooth(aes(x = month, y = offset_ratio), # global trend
color = "black", linetype = "dashed",
se = FALSE, span = 0.3) +
labs(
title = "Monthly blame Over Time by Party with Global Trend",
x = "Month",
y = "Blame pr sentence"
) +
theme_minimal()## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : span too small. fewer data values than degrees of freedom.## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : at 2022.4## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : radius 6.9444e-07## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : all data on boundary of neighborhood. make span bigger## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 2022.4## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 0.00083333## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 1## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : at 2022.6## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : radius 6.9444e-07## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : all data on boundary of neighborhood. make span bigger## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 6.9444e-07## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : zero-width neighborhood. make span bigger
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : zero-width neighborhood. make span bigger## Warning: Failed to fit group 2.
## Caused by error in `predLoess()`:
## ! NA/NaN/Inf in foreign function call (arg 5)## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
#Do analysis `
model_intercept_18 <- glmmTMB(total_blame ~ 1 + (1|party) + (1|month_index)+ offset(log(total_nr_sentences)),
ziformula = ~1,
family = nbinom1,
data = monthly_blame_18)
summary(model_intercept_18)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ 1 + (1 | party) + (1 | month_index) + offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame_18
##
## AIC BIC logLik -2*log(L) df.resid
## 3779.0 3800.2 -1884.5 3769.0 503
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.07798 0.2792
## month_index (Intercept) 0.05893 0.2428
## Number of obs: 508, groups: party, 13; month_index, 48
##
## Dispersion parameter for nbinom1 family (): 2.15
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.74911 0.08866 -31.01 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.45 2067.70 -0.01 0.992model_in_gov_18 <- glmmTMB(total_blame ~ in_gov + (1|party) + (1|month_index)+ offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = monthly_blame_18)
summary(model_in_gov_18)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + (1 | party) + (1 | month_index) + offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame_18
##
## AIC BIC logLik -2*log(L) df.resid
## 3761.5 3786.8 -1874.7 3749.5 502
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.05138 0.2267
## month_index (Intercept) 0.05870 0.2423
## Number of obs: 508, groups: party, 13; month_index, 48
##
## Dispersion parameter for nbinom1 family (): 2.05
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.72448 0.07605 -35.83 < 2e-16 ***
## in_govTRUE -0.35780 0.07945 -4.50 6.69e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.67 2307.87 -0.009 0.993model_block_18 <- glmmTMB(total_blame ~ in_gov + block + (1|party) + (1|month_index) + offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = monthly_blame_18)
summary(model_block_18)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + block + (1 | party) + (1 | month_index) +
## offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame_18
##
## AIC BIC logLik -2*log(L) df.resid
## 3758.6 3788.2 -1872.3 3744.6 501
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.03331 0.1825
## month_index (Intercept) 0.05845 0.2418
## Number of obs: 508, groups: party, 13; month_index, 48
##
## Dispersion parameter for nbinom1 family (): 2.05
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.87114 0.08823 -32.54 < 2e-16 ***
## in_govTRUE -0.35356 0.07766 -4.55 5.3e-06 ***
## blockRight 0.26927 0.11060 2.43 0.0149 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.98 2696.65 -0.008 0.993model_wing_18 <- glmmTMB(total_blame ~ in_gov + wing + (1|party) + (1|month_index)+ offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = monthly_blame_18)
summary(model_wing_18)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + wing + (1 | party) + (1 | month_index) +
## offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame_18
##
## AIC BIC logLik -2*log(L) df.resid
## 3759.8 3789.5 -1872.9 3745.8 501
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.03859 0.1964
## month_index (Intercept) 0.05869 0.2423
## Number of obs: 508, groups: party, 13; month_index, 48
##
## Dispersion parameter for nbinom1 family (): 2.05
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.9344 0.1237 -23.726 < 2e-16 ***
## in_govTRUE -0.3375 0.0791 -4.267 1.98e-05 ***
## wingWing 0.2757 0.1368 2.015 0.0439 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.61 2240.80 -0.01 0.992model_wing_block_18 <- glmmTMB(total_blame ~ in_gov + wing + block + (1|party) + (1|month_index) + offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = monthly_blame_18)
summary(model_wing_block_18)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + wing + block + (1 | party) + (1 | month_index) +
## offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame_18
##
## AIC BIC logLik -2*log(L) df.resid
## 3750.0 3783.8 -1867.0 3734.0 500
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.01110 0.1054
## month_index (Intercept) 0.05819 0.2412
## Number of obs: 508, groups: party, 13; month_index, 48
##
## Dispersion parameter for nbinom1 family (): 2.06
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.17631 0.09475 -33.52 < 2e-16 ***
## in_govTRUE -0.28184 0.07696 -3.66 0.00025 ***
## wingWing 0.34629 0.08270 4.19 2.82e-05 ***
## blockRight 0.33845 0.07221 4.69 2.77e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.57 2195.36 -0.01 0.992model_wing_block_int_18 <- glmmTMB(total_blame ~ in_gov + wing + block + wing*block + (1|party) + (1|month_index) + offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = monthly_blame_18)
summary(model_wing_block_int)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + wing + block + wing * block + (1 | party) +
## (1 | month_index) + offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 18414.0 18465.8 -9198.0 18396.0 2321
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.05221 0.2285
## month_index (Intercept) 0.05299 0.2302
## Number of obs: 2330, groups: party, 13; month_index, 274
##
## Dispersion parameter for nbinom1 family (): 3.14
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.81219 0.22947 -12.255 <2e-16 ***
## in_govTRUE -0.36992 0.01987 -18.618 <2e-16 ***
## wingWing 0.06676 0.25273 0.264 0.792
## blockRight -0.13827 0.28049 -0.493 0.622
## wingWing:blockRight 0.24231 0.32004 0.757 0.449
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.64 1105.78 -0.02 0.984model_wing_block_int_2_18 <- glmmTMB(total_blame ~ in_gov + block + in_gov*block + (1|party) + (1|month_index) + offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = monthly_blame_18)
summary(model_wing_block_int_2)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + block + in_gov * block + (1 | party) +
## (1 | month_index) + offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 18401.7 18447.7 -9192.8 18385.7 2322
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.06153 0.2481
## month_index (Intercept) 0.05271 0.2296
## Number of obs: 2330, groups: party, 13; month_index, 274
##
## Dispersion parameter for nbinom1 family (): 3.12
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.748347 0.105649 -26.014 < 2e-16 ***
## in_govTRUE -0.444119 0.028757 -15.444 < 2e-16 ***
## blockRight -0.001707 0.144255 -0.012 0.990559
## in_govTRUE:blockRight 0.177978 0.050112 3.552 0.000383 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.37 956.14 -0.022 0.982model_wing_block_int_3_18 <- glmmTMB(total_blame ~ in_gov + wing + block + in_gov*block + (1|party) + (1|month_index) + offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = monthly_blame_18)
summary(model_wing_block_int_3)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + wing + block + in_gov * block + (1 | party) +
## (1 | month_index) + offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 18401.7 18453.5 -9191.9 18383.7 2321
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.05280 0.2298
## month_index (Intercept) 0.05272 0.2296
## Number of obs: 2330, groups: party, 13; month_index, 274
##
## Dispersion parameter for nbinom1 family (): 3.12
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.93526 0.16169 -18.154 < 2e-16 ***
## in_govTRUE -0.44354 0.02876 -15.423 < 2e-16 ***
## wingWing 0.22576 0.15634 1.444 0.148715
## blockRight 0.02874 0.13591 0.211 0.832531
## in_govTRUE:blockRight 0.17926 0.05013 3.576 0.000349 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.54 1041.90 -0.021 0.984model_wing_block_int_4_18 <- glmmTMB(total_blame ~ in_gov + wing + block + wing*block + in_gov*wing + in_gov*block + (1|party) + (1|month_index) + offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = monthly_blame_18)## dropping columns from rank-deficient conditional model: in_govTRUE:wingWingsummary(model_wing_block_int_4)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + wing + block + wing * block + in_gov *
## wing + in_gov * block + (1 | party) + (1 | month_index) +
## offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 18404.6 18467.9 -9191.3 18382.6 2319
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.04898 0.2213
## month_index (Intercept) 0.05271 0.2296
## Number of obs: 2330, groups: party, 13; month_index, 274
##
## Dispersion parameter for nbinom1 family (): 3.12
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.78716 0.22250 -12.526 < 2e-16 ***
## in_govTRUE -0.43716 0.03276 -13.344 < 2e-16 ***
## wingWing 0.04535 0.24511 0.185 0.853202
## blockRight -0.19506 0.27226 -0.716 0.473716
## wingWing:blockRight 0.29039 0.31057 0.935 0.349779
## in_govTRUE:wingWing -0.02709 0.05868 -0.462 0.644287
## in_govTRUE:blockRight 0.17416 0.05273 3.303 0.000956 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.55 1047.29 -0.021 0.984anova(model_intercept_18, model_in_gov_18, model_wing_18, model_wing_block_18, model_wing_block_int_18,model_wing_block_int_2_18,model_wing_block_int_3_18,model_wing_block_int_4_18)## Data: monthly_blame_18
## Models:
## model_intercept_18: total_blame ~ 1 + (1 | party) + (1 | month_index) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_in_gov_18: total_blame ~ in_gov + (1 | party) + (1 | month_index) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_wing_18: total_blame ~ in_gov + wing + (1 | party) + (1 | month_index) + , zi=~1, disp=~1
## model_wing_18: offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_wing_block_18: total_blame ~ in_gov + wing + block + (1 | party) + (1 | month_index) + , zi=~1, disp=~1
## model_wing_block_18: offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_wing_block_int_2_18: total_blame ~ in_gov + block + in_gov * block + (1 | party) + , zi=~1, disp=~1
## model_wing_block_int_2_18: (1 | month_index) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_wing_block_int_18: total_blame ~ in_gov + wing + block + wing * block + (1 | party) + , zi=~1, disp=~1
## model_wing_block_int_18: (1 | month_index) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_wing_block_int_3_18: total_blame ~ in_gov + wing + block + in_gov * block + (1 | party) + , zi=~1, disp=~1
## model_wing_block_int_3_18: (1 | month_index) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_wing_block_int_4_18: total_blame ~ in_gov + wing + block + wing * block + in_gov * , zi=~1, disp=~1
## model_wing_block_int_4_18: wing + in_gov * block + (1 | party) + (1 | month_index) + , zi=~1, disp=~1
## model_wing_block_int_4_18: offset(log(total_nr_sentences)), zi=~1, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df
## model_intercept_18 5 3779.0 3800.2 -1884.5 3769.0
## model_in_gov_18 6 3761.5 3786.8 -1874.7 3749.5 19.5800 1
## model_wing_18 7 3759.8 3789.5 -1872.9 3745.8 3.6179 1
## model_wing_block_18 8 3750.0 3783.8 -1867.0 3734.0 11.8497 1
## model_wing_block_int_2_18 8 3759.9 3793.7 -1871.9 3743.9 0.0000 0
## model_wing_block_int_18 9 3750.3 3788.4 -1866.1 3732.3 11.5863 1
## model_wing_block_int_3_18 9 3752.0 3790.1 -1867.0 3734.0 0.0000 0
## model_wing_block_int_4_18 10 3751.9 3794.2 -1865.9 3731.9 2.1233 1
## Pr(>Chisq)
## model_intercept_18
## model_in_gov_18 9.647e-06 ***
## model_wing_18 0.0571627 .
## model_wing_block_18 0.0005767 ***
## model_wing_block_int_2_18 1.0000000
## model_wing_block_int_18 0.0006644 ***
## model_wing_block_int_3_18 1.0000000
## model_wing_block_int_4_18 0.1450711
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(model_in_gov_18, model_wing_block_18, model_wing_block_int_18)## Data: monthly_blame_18
## Models:
## model_in_gov_18: total_blame ~ in_gov + (1 | party) + (1 | month_index) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_wing_block_18: total_blame ~ in_gov + wing + block + (1 | party) + (1 | month_index) + , zi=~1, disp=~1
## model_wing_block_18: offset(log(total_nr_sentences)), zi=~1, disp=~1
## model_wing_block_int_18: total_blame ~ in_gov + wing + block + wing * block + (1 | party) + , zi=~1, disp=~1
## model_wing_block_int_18: (1 | month_index) + offset(log(total_nr_sentences)), zi=~1, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df
## model_in_gov_18 6 3761.5 3786.8 -1874.7 3749.5
## model_wing_block_18 8 3750.0 3783.8 -1867.0 3734.0 15.4675 2
## model_wing_block_int_18 9 3750.3 3788.4 -1866.1 3732.3 1.7052 1
## Pr(>Chisq)
## model_in_gov_18
## model_wing_block_18 0.0004378 ***
## model_wing_block_int_18 0.1916153
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- model_wing_block_18 is the better model for predicting amount of blame within recent years
#check assumption of best model
model_check(model_wing_block_18)
## chisq ratio rdf p
## 458.0709762 0.9161420 500.0000000 0.9104718#Plot best model
summary(model_wing_block_18)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ in_gov + wing + block + (1 | party) + (1 | month_index) +
## offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame_18
##
## AIC BIC logLik -2*log(L) df.resid
## 3750.0 3783.8 -1867.0 3734.0 500
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.01110 0.1054
## month_index (Intercept) 0.05819 0.2412
## Number of obs: 508, groups: party, 13; month_index, 48
##
## Dispersion parameter for nbinom1 family (): 2.06
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.17631 0.09475 -33.52 < 2e-16 ***
## in_govTRUE -0.28184 0.07696 -3.66 0.00025 ***
## wingWing 0.34629 0.08270 4.19 2.82e-05 ***
## blockRight 0.33845 0.07221 4.69 2.77e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.57 2195.36 -0.01 0.992sjPlot::plot_model(model_wing_block_18,
type = "pred",
terms = c("in_gov", "wing"),
condition = c(total_nr_sentences = 100),
axis.title = c("In Government", "Predicted Blame Counts pr 100 Uttered Sentences"),
title = "Predicted Blame by In Government and Wing in period 2019-2022",
ci.lvl = NA,
axis.lim = c(3,6.5)
) +
theme_minimal()## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
sjPlot::plot_model(model_wing_block_18,
type = "pred",
terms = c("block"),
condition = c(total_nr_sentences = 100),
axis.title = c("Block", "Predicted Blame Counts pr 100 Uttered Sentences"),
title = "Predicted Blame by Block in period 2019-2022",
ci.lvl = NA,
axis.lim = c(3,6.5)
) +
theme_minimal()## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
While the general amount of blame has not changed according to this paper, it seems like the polarization has become more clear than before with also wing and block in addition to in_gov being significant predictors of
(Intercept) -3.16130 0.09924 -31.86 < 2e-16 in_govTRUE -0.27310 0.08655 -3.16 0.0016 wingWing 0.34538 0.08697 3.97 7.15e-05 blockRight 0.32434 0.07422 4.37 1.24e-05 *
fixef(model_wing_block_int_18)##
## Conditional model:
## (Intercept) in_govTRUE wingWing
## -3.0513 -0.3243 0.1955
## blockRight wingWing:blockRight
## 0.1691 0.2199
##
## Zero-inflation model:
## (Intercept)
## -21.9#time as predictor of blame in recent years
model_blame_time_18 <- glmmTMB(total_blame ~ month_index + (1|party) + (1|in_gov) + offset(log(total_nr_sentences)),
ziformula = ~1,
family = nbinom2,
data = monthly_blame_18)
summary(model_blame_time_18)## Family: nbinom2 ( log )
## Formula:
## total_blame ~ month_index + (1 | party) + (1 | in_gov) + offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: monthly_blame_18
##
## AIC BIC logLik -2*log(L) df.resid
## 3927.2 3952.6 -1957.6 3915.2 502
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.05319 0.2306
## in_gov (Intercept) 0.03281 0.1811
## Number of obs: 508, groups: party, 13; in_gov, 2
##
## Dispersion parameter for nbinom2 family (): 7.85
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.036470 0.393107 -10.268 < 2e-16 ***
## month_index 0.004711 0.001420 3.318 0.000907 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.6 2430.8 -0.009 0.993#Plot this model
sjPlot::plot_model(model_blame_time_18,
type = "pred",
terms = c("month_index"),
condition = c(total_nr_sentences = 100),
axis.title = c("Month since Jan 2000", "Predicted Blame Counts pr 100 Uttered Sentences"),
title = "Evolution of the Use of Blame in period 2019-2022",
ci.lvl = NA,
axis.lim = c(3,6.5)
) +
theme_minimal()## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
#Now the trend is turning and in the last few years it has been going up
# Define election months
election_months <- as.yearmon(c("Feb 2005", "Nov 2007", "Sep 2011", "Jun 2015", "Jun 2019"))
monthly_blame$time_to_election <- ifelse(monthly_blame$month %in% election_months, 0, NA )
# Assume monthly_blame$month is also yearmon
# If not, convert:
# monthly_blame$month <- as.yearmon(monthly_blame$month)
# Function to compute time to nearest election
time_to_nearest_election <- function(month, elections, window = 12) {
# compute differences in months
diffs <- as.numeric((month - elections) * 12)
# find the nearest election
closest <- diffs[which.min(abs(diffs))]
# return NA if farther than ±window months
if (abs(closest) > window) {
return(NA_real_)
} else {
return(+closest) # negative before, positive after (flip sign for intuition)
}
}
# Apply function to each month
monthly_blame <- monthly_blame %>%
mutate(time_since_election = sapply(month, time_to_nearest_election, elections = election_months))election_monthly_blame <- monthly_blame %>%
filter(!is.na(time_since_election))election_model <- glmmTMB(total_blame ~ time_since_election + (1|party) + (1|in_gov) + offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = election_monthly_blame)
summary(election_model)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ time_since_election + (1 | party) + (1 | in_gov) +
## offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: election_monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 8256.9 8286.4 -4122.4 8244.9 1013
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.04144 0.2036
## in_gov (Intercept) 0.01971 0.1404
## Number of obs: 1019, groups: party, 11; in_gov, 2
##
## Dispersion parameter for nbinom1 family (): 5.64
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.902869 0.119574 -24.277 <2e-16 ***
## time_since_election 0.003809 0.001570 2.426 0.0153 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -20.66 1416.04 -0.015 0.988election_model_quad <- glmmTMB(total_blame ~ scale(time_since_election) + I(scale(time_since_election)^2) + (1|party) + (1|in_gov) + offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = election_monthly_blame)
summary(election_model_quad)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ scale(time_since_election) + I(scale(time_since_election)^2) +
## (1 | party) + (1 | in_gov) + offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: election_monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 8253.7 8288.2 -4119.9 8239.7 1012
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.04131 0.2033
## in_gov (Intercept) 0.02029 0.1425
## Number of obs: 1019, groups: party, 11; in_gov, 2
##
## Dispersion parameter for nbinom1 family (): 5.58
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.93563 0.12162 -24.139 <2e-16 ***
## scale(time_since_election) 0.02767 0.01107 2.499 0.0124 *
## I(scale(time_since_election)^2) 0.02833 0.01246 2.274 0.0230 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -20.25 1510.77 -0.013 0.989election_model_triq <- glmmTMB(total_blame ~ scale(time_since_election) + I(scale(time_since_election)^2)+ + I(scale(time_since_election)^3) + (1|party) + (1|in_gov) + offset(log(total_nr_sentences)),
ziformula = ~1, # models extra zeros
family = nbinom1,
data = election_monthly_blame)
summary(election_model_triq)## Family: nbinom1 ( log )
## Formula:
## total_blame ~ scale(time_since_election) + I(scale(time_since_election)^2) +
## +I(scale(time_since_election)^3) + (1 | party) + (1 | in_gov) +
## offset(log(total_nr_sentences))
## Zero inflation: ~1
## Data: election_monthly_blame
##
## AIC BIC logLik -2*log(L) df.resid
## 8255.6 8295.0 -4119.8 8239.6 1011
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## party (Intercept) 0.04134 0.2033
## in_gov (Intercept) 0.02027 0.1424
## Number of obs: 1019, groups: party, 11; in_gov, 2
##
## Dispersion parameter for nbinom1 family (): 5.58
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.935307 0.121579 -24.143 <2e-16 ***
## scale(time_since_election) 0.017379 0.030430 0.571 0.5679
## I(scale(time_since_election)^2) 0.028209 0.012464 2.263 0.0236 *
## I(scale(time_since_election)^3) 0.005437 0.014982 0.363 0.7167
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -20.42 1539.11 -0.013 0.989ggplot(election_monthly_blame, aes(x = time_since_election, y=offset_ratio, color = party)) +
geom_smooth()+
geom_smooth(aes(x = time_since_election, y=offset_ratio), color = "black", linetype = "dashed")## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
model_check(election_model_quad)
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
## chisq ratio rdf p
## 996.1515889 0.9843395 1012.0000000 0.6328048#plot best model
sjPlot::plot_model(election_model_quad,
type = "pred",
terms = "time_since_election[all]",
condition = c(total_nr_sentences = 100),
axis.title = c("Months Since Nearest Election", "Predicted Blame Counts pr 100 Uttered Sentences"),
title = "Predicted Counts of blame by time since nearest election ",
ci.lvl = NA,
show.legend = TRUE
) +
theme_minimal()## Model uses a transformed offset term. Predictions may not be correct.
## It is recommended to fix the offset term using the `condition` argument,
## e.g. `condition = c(party = 1)`.
## You could also transform the offset variable before fitting the model
## and use `offset(party)` in the model formula.